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Abstract 

Direct  numerical  simulations  of  compressible  fluid  flow  are  performed  for 
subsonic  and  supersonic  channel  flow  with  two  symmetrically  backward  fac- 
ing steps  and  for  a supersonic  compression  ramp  flow  field.  Spatial  deriva- 
tives are  represented  by  a central  scheme  of  high  order  difference  operators 
(N=2,4,6,8,...)  that  is  used  together  with  artificial  dissipation  of  order  N. 
A two-step  Richtmyer  scheme  is  employed  for  time  integration.  In  regions 
with  steep  gradients  flux-corrected  transport  (FCT)  according  to  Boris  and 
Book  is  applied.  Preliminary  results  are  presented  for  Mach  number  1.5  in 
case  of  the  channel  flow  and  further  results  for  Mach  number  2.84  in  case 
of  the  ramp  with  a ramp  angle  of  24  degrees. 

1.  Introduction 

The  Reynolds  Averaged  Navier-Stokes  equations  (RANS)  are  often  used  in 
connection  with  various  turbulence  models  to  model  technical  problems  of 
compressible  fluid  flow,  e.g.  [1].  But  the  aptitude  of  a particular  turbulence 
model  for  the  problem  to  solve  is  generally  not  known  beforehand  and  may 
even  be  unsatisfactory.  On  the  other  hand,  the  rapidly  growing  computer 
resources  offer  a promising  future  to  physically  more  realistic  mathematical 
models  like  Large-Eddy  Simulation  (LES)  and  Direct  Numerical  Simulation 
(DNS),  see  for  example  [2] [3].  Both  methods  require  high  order  discretiza- 
tion with  special  treatment  of  shocks  or  steep  gradients.  While  LES  still 
needs  some  modeling  to  account  for  the  spatially  underresolved  stresses, 
DNS  should  be  apt  to  resolve  all  turbulent  scales  without  the  support  of 
any  empirical  turbulence  modeling.  While  the  method  in  [3]  employs  al- 
ternating upwinding  using  compact  differences,  the  method  applied  in  this 


479 


480 


I.  KLUTCHNIKOV  AND  J.  BALLMANN 


paper  is  based  on  central  high  order  numerical  operators  for  interpola- 
tion and  approximation  in  space  in  conjunction  with  a two-step  Richtmyer 
scheme.  The  method  represents  an  extension  of  the  concept  proposed  in 
[4]  for  the  Burgers’  equation  and  the  Euler  equations  to  the  system  of 
the  Navier-Stokes  equations  for  heat  conducting  compressible  fluid  flow  [5]. 
Flux-Corrected  Transport  (FCT)  according  to  Boris  and  Book  has  been 
implemented  to  account  for  steep  gradients,  e.g.  in  the  presence  of  shocks. 

2.  Numerical  Method 

2.1.  INTERIOR  POINTS  OF  THE  SOLUTION  DOMAIN 

Denoting  by  U the  set  of  conserved  variables  p,  pu,  pv,  pw , petotai  and  by 
Fc  the  fluxes  corresponding  to  the  inviscid  and  by  Fd  the  dissipative  part  of 
the  fluxes  and  leaving  out  body  forces  and  body  energy  supply  the  system 
of  Navier-Stokes  equations  reads 


Ut  + j^F^j^Ff  (1) 

r= 1 t— 1 

The  equations  are  integrated  with  respect  to  time  employing  a second  or- 
der Richtmyer  scheme.  Space  discretization  is  performed  with  central  high 
order  numerical  operators  [4]  [5].  The  numerical  scheme  requires  artificial 
dissipation  of  the  highest  order  of  the  scheme  for  numerical  stability  in  the 
sense  of  von  Neumann  [5].  Using  lower  index  ir  for  discretization  in  the 
space  directions  r,  r = 1,2,3  and  upper  index  n for  time  stepping  the  set 
of  discretized  equations  reads 

= LrU?T+1/2  - 0.5A rAX+1/2RrK+l/2 

3 

jrrn  \ ' \ r t viCn-\-l/2  j-,crc-tT/2\ 

~ U»  r+l/2  *rir-l/2> 

+ (C+l/2  - *£-,/ 2>  - (Sr^+1/2  - SrUl_in)} 

ir  = {i,j,k}  , At/Axr=:Xr  , r = l,2, 3 (2) 

Therein  Lr , Rr)  Sr  represent  the  aforementioned  high  order  operators,  Lr 
serving  for  interpolation,  Rr  for  approximation  and  Sr  for  artificial  dissi- 
pation. Ar  is  the  Jacobian.  With  coefficients  am  in  Lr,  bm  in  Rr , and  drn 
in  Sr  the  operators  of  order  N are  defined  as  follows, 

N/2 

LrUir+ 1/2  = £ ( — 1 )m+1am(Uir+m  + Uir-m+l), 

m—  1 
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(3) 
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N/2 

RrUir+ 1/2  = E ("I r+1bm(Uir+m  ~ Uir-m+ l),  (4) 

m=l 

N/  2 

SrUir+\/2  = Q:ir  + l/2  EM)  + dm(Uir+m  + [Zjr-TO+l))  (5) 

m=l 

with  air+i/2  a prescribed  factor  for  each  cell  index.  The  expression  for  the 
artificial  dissipation  of  order  N is 

f)Nrj 

-a-^Ax?  = -(SrUir+l/2  - SrUir_l/2) Ax"1,  (6) 

with  a a number  to  be  chosen.  Optionally,  in  domains  with  steep  gradients 
FCT  is  applied  in  a self-controlled  fashion  with  the  following  expressions 
for  the  anti-diffusive  flux  (index  ad)  and  the  corrected  flux  (index  cor) 


pc(ad)  _ p,c(/i)n+l/2  _ pC.(l)n+ 1/2 
PTir  + 1/2  - rrir+ 1/2  PriT  + 1/2  ’ 


pc(cor)n+l/2 

Pri,+l/2 


0+1/2 


p,c(arf) 

'rrir+l/2' 


(7) 


The  coefficients  am,bm  and  dm  in  Eq.  (3),  (4),  and  (5)  are  depending  on 
the  order  N which  is  chosen  for  the  solution  [4]  [5]. 


2.2.  BOUNDARY  POINTS 

At  solid  walls  the  no  slip  condition  is  prescribed  and  walls  are  assumed 
thermally  adiabatical.  The  central  scheme  of  order  N makes  use  of  N/2 
fictitious  points  by  mirror  principle.  For  solution  points  at  artificial  bound- 
aries marking  the  boundary  of  the  computational  domain,  e.g.  at  inflow  and 
outflow,  also  N/2  fictitious  points  are  needed.  Different  conditions  are  to  be 
distinguished  in  the  fictitious  points  there  for  subsonic  and  supersonic  flow. 
For  subsonic  inflow  and  outflow  part  of  the  values  in  the  fictitious  points 
are  set  using  Riemann  invariants  according  to  a concept  of  local  simple 
waves  in  the  sense  of  gasdynamics  while  the  other  part  of  the  values  are 
extrapolated.  This  way  non-physical  reflections  from  the  artificial  boundary 
are  mostly  suppressed. 


2.3.  APPLICATION  PROCEDURE  AND  VALIDATION 

The  examples  we  discuss  in  this  paper  would  represent  two-dimensional 
flow  fields  in  case  of  laminar  flow.  The  calculation  starts  as  for  laminar 
flow.  On  the  inflow  boundary  at  a wall  a starting  boundary  layer  with 
parabolic  velocity  profile  is  prescribed  which  may  become  later  on  changed 
as  part  of  the  solution  because  of  locally  subsonic  inflow.  At  walls  the 
correct  boundary  conditions  are  introduced  from  the  first  iteration  step. 
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First,  some  thousands  of  time  steps  of  the  solution  are  performed  for  two- 
dimensional  flow.  Then  during  a certain  number  of  time  steps  a disturbance 
consisting  of  two  additional  velocity  components  in  the  cross-plane  of  the 
main  flow  direction  is  superposed  to  the  boundary  condition  at  the  entrance 
of  the  solution  domain  which  makes  the  flow  three-dimensional.  Thereafter 
this  disturbance  is  removed. 

Numerous  and  extended  validation  tests  which  are  not  presented  in  the  pa- 
per have  been  conducted  for  the  method  [5]  [11],  e.g.  comparison  was  made 
with  analytical  solutions  of  the  Burgers’  equation,  of  the  ID  Riemann  prob- 
lem for  the  Euler  equations,  and  of  the  mean  velocity  distributions  in  the 
viscous  layer  and  the  logarithmic  layer  from  turbulent  boundary  layer  the- 
ory. Results  for  compressible  flow  at  Mach  number  M=0.2  in  a straight 
channel  at  Reynolds  number  Re=1750  have  been  compared  with  DNS  re- 
sults [6]  [7]  and  experimental  results  at  comparable  Reynolds  numbers  [8]  [9]. 
Numerical  results  for  a Mach  0.2  flow  with  Re=2600  through  a channel  with 
a backward  facing  step  were  compared  with  experimentally  determined  val- 
ues [10]  including  the  Reynolds  stresses  at  different  positions  behind  the 
step  [5]  [11].  A kind  of  self- validation  is  the  check  of  asymptotic  convergence 
with  increasing  the  order  N first  and  then  refining  the  grid.  This  has  been 
done  e.g.  in  [11]  for  the  flow  through  a channel  with  a straight  axis  and  a 
symmetric  jump  of  the  channel  height  from  h on  the  upstream  side  to  3h  on 
the  downstream  side.  Data  for  Mach  and  Reynolds  numbers  were  M=0.6 
and  Re=103 4.  Asymptotic  convergence  was  found  for  N=8. 

3.  Simulations  of  Supersonic  Flow  Fields 

A Mach  1.5  flow  through  the  aforementioned  channel  with  two  symmetrical 
backwards  facing  steps  has  been  simulated.  Static  pressure  in  the  flow  just 
before  and  just  behind  the  steps  has  the  same  value  p=l.  The  Reynolds 
number  based  on  the  step  height  h is  Re=104.  The  length  of  the  solution 
domain  from  the  steps  to  the  outflow  end  is  16h.  Number  of  grid  cells  in 
this  domain  is  256x61x41.  Flow  is  from  left  to  right.  Computation  is  started 
with  given  supersonic  inflow  upstream  of  the  steps  including  a boundary 
layer  of  thickness  0.15h  at  the  upper  and  lower  wall  with  a parabolic  veloc- 
ity profile.  In  Fig.(l)  results  of  the  density  distribution  are  shown  for  orders 
from  N=2  up  to  N=16  for  the  same  physical  time  t=20  which  means  20000 
time  steps.  The  flow  field  is  not  yet  fully  established.  Comparing  the  contour 
lines  for  the  different  orders  one  recognizes  that  contrary  to  the  subsonic 
M=0.6  case  the  order  N=8  seems  not  sufficient  for  asymptotic  convergence. 
Before  the  order  step  from  N=14  to  16  the  solution  exhibits  decreasing  but 
still  remarkable  changes  with  increasing  order.  The  long  time  simulation 
was  then  conducted  with  order  16.  In  Fig. (2)  some  results  obtained  so  far 
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are  presented.  The  instantaneous  contour  lines  of  Mach  number,  pressure 
and  density  distributions  for  t =70  show  on  the  inflow  side  the  typical  jet 
behavior  for  the  case  of  adapted  pressure.  Then  the  jet  becomes  declining 
more  and  more,  exhibiting  vortical  structures,  local  supersonic  pockets  at 
vortices  and  shocklets.  The  ranges  of  Mach  number  are  from  M=0.1  to 
M=2.8  and  of  pressure  and  density  from  p=0.3  to  p=2.3  and  p=0.4  to  1.7, 
respectively.  Pressure  fluctuations  for  the  time  interval  from  t=40  to  t=100 
have  been  recorded  in  points  P : x\  = 3 A,  aq  = 2 h,  = 0 and  P : x\  = 15 h, 
X2  = 2 h,  ,x'3  = 0.  The  graphs  and  frequency  spectra  do  not  yet  represent 
fully  developed  turbulent  flow.  It  seems  that  the  computational  domain  is 
not  long  enough  to  observe  the  full  transition  region.  Further  investiga- 
tions are  necessary.  As  a second  example  the  supersonic  flow  over  a 24° 
ramp  is  revisited  [11]  Mach  number  and  Reynolds  number  are  M=2.84  and 
Re=107.  According  to  experiments  by  Settles  [12]  [13] . results  are  depicted 
in  Fig.  (3).  The  length  of  the  separation  zone  is  reproduced  fairly  well  by 
the  simulation,  and  the  three-dimensional  character  of  the  flow  is  obvious. 
The  pressure  signature  in  the  middle  of  the  corner  line  is  within  a range 
of  p=1.9  to  p=3.3.  The  mean  pressure  is  overpredicted  in  the  separation 
zone  but  fits  fairly  well  the  experimental  values  outside.  The  skin  friction 
coefficient  fits  best  the  earlier  experimental  values  found  by  Settles  [12]  . 
That  agrees  with  recent  results  [2]. 


References 

1.  Rodi  W.,  Mansour  N.N.  (1993)  Low  Reynolds  number  k-e  modeling  with  the  aid  of 
direct  simulation  data,  J.  Fluid  Mech.  250,  pp.  509-529 

2.  Rizzetta,  D.P.  and  Visbal,  M.R.  (2001)  Large-Eddy  Simulation  of  Supersonic 
Compression-Ramp  Flows,  AIAA  paper  2001-2858 

3.  Fezer  A.,  Kloker  M.  (1998)  Transition  Processes  in  Mach  6.8  Boundary  Layers  at 
Varying  Temperature  Conditions  Investigated  by  Spatial  Direct  Numerical  Simula- 
tion, Notes  on  Numerical  Fluid  Mechanics,  Vol.  no.  72,  pp.  138-145 

4.  Zalesak  S.T.  (1984)  A Physical  Interpretation  of  the  Richtmyer  Two-Step  Lax- 
Wendroff  Scheme,  and  its  Generalization  to  Higher  Spatial  Order” , Adv.  In  Comp. 
Meth.  For  Part.  Diff.  Eq.,  Publ.  IMACS,  pp.  491-496 

5.  Klutchnikov,  I.  (1998)  Direct  numerical  simulation  of  turbulent  compressible  fluid 
flow.  Habilitation  thesis  (in  Russian),  Russian  Academy  of  Sciences,  Moscow 

6.  Kim  J.,  Moin  P.,  Moser  R.D.  (1987)  Turbulence  Statistics  in  Fully  Developed  Chan- 
nel Flow  at  Low  Reynolds  Number,  J.  of  Fluid  Mech.  177,  pp.  133-166 

7.  Spalart  P.R.  (1988)  Direct  simulation  of  a turbulent  boundary  layer  up  to  R=1410, 
J.  of  Fluid  Mech.  187,  pp.  61-98 

8.  Whan  G.A.,  Rothfus  R.R.  (1959)  Amer.  Inst.  Chem.  Eng.  J.,  Vol.  no.  5,  N2, 
pp.  205-208 

9.  Patel  V.C.,  Head  M.R.  (1969)  Some  Observations  on  Skin  Friction  and  Velocity 
Profiles  in  Fully  Developed  Pipe  and  Channel  Flows,  J.  Fluid  Mech.,  Vol.  no.  38, 
Nl,  pp.  181-201 

10.  Komarov  P.L.,  Polyakov  A.F.  (1996)  Investigation  About  the  Characteristics  of  Tur- 
bulence and  Heat  Exchange  in  a Channel  Behind  a Backward  Facing  Step.  Preprint 
IVTAN  N2-396,  Moscow,  70p.  (in  Russian) 


484 


I.  KLUTCHNIKOV  AND  J.  BALLMANN 


11.  Klutchnikov  I.,  Ballmann  J.  (2001)  DNS  of  Turbulent  Compresible  Fluid  Flow  with 
a High  Order  Difference  Method,  AIAA  paper  2001-2542 

12.  Settles  G.S.,  Vas  I.E.,  BogdonofF  S.M.  (1976)  Details  of  a Shock-Separated  Turbu- 
lent Boundary  Layer  at  a Compression  Corner,  AIAA  Journal,  Vol.  14,  No.  12, 
pp.  1709-1715 

13.  Settles  G.S.,  Dodson  L.J.  (1991)  Hypersonic  Shock/Boundary-Layer  Interaction 
Database,  AIAA  paper  1991-1763 
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Figure  1.  Supersonic  flow  (M=1.5,  Ee=104)  through  a channel  with  two  backwards  fac- 
ing steps.  Asymptotic  convergence  check  of  the  instantaneous  density  field  with  increasing 
order,  t=20 
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Figure  2.  Supersonic  channel  flow  (M=1.5,  Re=104)  , instantaneous  contours  of  Mach 
number,  pressure  and  density  for  time  t=70  and  signatures  of  pressure  over  time  from 
t=40  to  t=100  in  points  Pi:  X\  = 3 h,  xi  = 2 h,  xz  = 0 and  P2:  x\  = 15 h,  xi  — 2 h, 
xs  = 0.  Lower  figures:  respective  Fourier  spectra 
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Figure  3.  Supersonic  flow  along  a 24°  compression  ramp  (M=2.84,  Re=107),  isosurfaces 
of  M=2.74.  2.5  and  2.05  and  p=2.58  and  1.54  (upper  figures),  graph  of  pressure  signature 
in  the  middle  point  of  the  corner  line  for  the  time  interval  from  t=30  to  t=60  and  the 
respective  Fourier  spectrum  (middle  figures)  and  mean  wall  distributions  of  pressure  p 
and  skin  friction  coefficient  ct  (lower  figures).  Symbols  represent  experimental  values 


